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Abstract 

We present results of numerical simulations of the 2+ld Nambu - Jona-Lasinio model 
with a non-zero baryon chemical potential /i including the effects of a diquark source 
term. Diquark condensates, susceptibilities and masses are measured as functions 
of source strength j. The results suggest that diquark condensation does not take 
place in the high density phase \x > /i c , but rather that the condensate scales non- 
analytically with j implying a line of critical points and long range phase coherence. 
Analogies are drawn with the low temperature phase of the 2d XY model. The 
spectrum of the spin-| sector is also studied yielding the quasiparticle dispersion 
relation. There is no evidence for a non-zero gap; rather the results are characteristic 
of a normal Fermi liquid with Fermi velocity less than that of light. We conclude that 
the high density phase of the model describes a relativistic gapless thin film BCS 
superfluid. 

PACS: ll.10.Kk, 11.30.Fs, 11.15.Ha, 21.65.+f, 67.70.+n 

Keywords: Monte Carlo simulation, chemical potential, diquark condensate, super- 
fluidity 



1 



1 Introduction 



Spontaneous symmetry breaking in particle physics was pre-dated by the Bardeen- 
Cooper-Schrieffer (BCS) mechanism for superconductivity in metals at low tempera- 
ture 0, which predicts a ground state in which a macroscopic fraction of the electrons 
in the vicinity of the Fermi surface reside in spin-0 bound states known as Cooper 
pairs. In field theoretic terms the Cooper pairs form a condensate which alter the 
symmetry of the ground state; in the case of superconductivity U(l) electromagnetic 
gauge invariance is spontaneously broken, leading to the Meissner effect, ie. exclusion 
of magnetic field from a superconducting sample due to surface screening currents. 
The ideas of BCS have been incorporated into particle physics in two distinct direc- 
tions. Firstly, particle - anti-particle pair condensation (ipip) ^ was suggested as a 
means of breaking the global chiral symmetries responsible for keeping fermion masses 
small; Goldstone's theorem then predicts light weakly-interacting bosonic states which 
can be identified with pions, whose masses are considerably less than the nucleon, the 
lightest strongly-interacting fermion. The resulting model provides a reasonable de- 
scription of low-energy strong- interaction phenomena [[|. Secondly, condensation of 
an elementary Higgs field has of course been invoked as a mechanism for electroweak 
gauge symmetry breaking, imbuing gauge bosons (as well as fermionic matter fields) 
with non-zero mass in precise analogy with the Meissner effect. 

In recent years the BCS mechanism has returned to particle physics in a new 
guise in the context of the fundamental theory of the strong interaction, QCD, at 
high density. For baryon charge densities % ~ 0(l)fm -3 , it is believed that chiral 
symmetry is restored and nucleons dissociate into quarks. The resulting ground 
state is thought to be quark matter in which to first approximation the dominant 
degrees of freedom are relativistic degenerate quarks forming a Fermi sphere with 
Fermi momentum kp ~ 350 — 400MeV. Such conditions are conceivably found in the 
cores of neutron stars. However, since the force between quarks due to eg. one-gluon 
exchange is attractive, this simple picture is unstable with respect to a BCS scenario 
in which condensation of diquark pairs occurs Q . Since the qq wavefunction is gauge 
non-singlet, the resulting ground state renders some or all of the gluons massive, an 
effect known as "color superconductivity" . The resulting dynamically generated mass 
scale or "gap" A is predicted to be O(100)MeV |J, and hence comparable with the 
consituent quark scale. 

Unfortunately, theoretical studies of color superconductivity are to date limited 
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to perturbative and self-consistent methods 0, |, |J; there is no systematic method 
of performing non-perturbative QCD calculations in the high density regime because 
of the notorious "sign problem", ie. the measure of the Euclidean path integral 
becomes complex once baryon chemical potential /i ^ 0, making the importance 
sampling techniques traditionally used in numerical simulations of lattice gauge theory 
ineffective. There are, however, strongly-interacting model field theories where this 
difficulty can be circumvented. One is QCD with just two colors, in which qq baryons 
and qq mesons fall into multiplets related by enhanced global symmetries. Some of 
the multiplets contain Goldstone bosons, so that the methods of chiral perturbation 
theory can be applied |7|. Baryonic matter ub > forms for chemical potential 
^ [0, H, [|. The resulting ground state is a superfluid of light but strongly 

bound qq states which form a Bose-Einstein condensate. There is no Fermi surface 
in this regime, and the BCS description is inappropriate. 

Another possibility, to be studied in the present paper, are four-fermi models such 
as the Nambu - Jona-Lasinio (NJL) model ||, in which it can be shown that the 
effects of adding "conjugate quarks" q c to make the path integral real and positive 
have little impact on the physics: at low density light Goldstone states arising as 
a result of chiral symmetry breaking can only form in mesonic qq channels, whereas 



baryonic qq c bound states remain massive, ie. at the constituent quark scale |T0[ . This 
means that unlike in physical three-color QCD with conjugate quarks, it is possible 
for simulations of four-fermi models to maintain a separation of scales between m v 
and the critical /x c at which chiral symmetry is restored and baryonic matter induced 
into the ground state Jill. ^ n this case the model does not reproduce any of the 
physics of confinement, but has a Fermi surface for tib > 0. 

Because the qq interaction is attractive, diquark condensation is expected in the 
high density phase of each of the models described above. In both cases, however, 
the relevant (qq) ^ is gauge singlet, meaning that the ground state is not supercon- 
ducting, but rather superfluid. In field theoretic terms, a superfluid forming by BCS 
condensation is characterised by a ground state which does not respect a global sym- 
metry of the underlying action, in this case the U(l) corresponding to baryon number, 
which is thus no longer a good quantum number. Fermionic excitations above the 
Fermi surface are a superposition of particle and hole states, and require energy > A 
to excite. Finally, because a continuous global symmetry is spontaneously broken, 
Goldstone's theorem applies and massless diquark states are expected in the exci- 
tation spectrum. Physically these result both in a long-ranged interaction between 
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the vortex excitations found in a rotating superfluid, and in propagating waves of 
temperature variation known as second sound |13|1 . 

The two known superfluids are liquid 4 He at kelvin and liquid 3 He at milli-kelvin 
temperatures. 4 He is a boson and is naturally treated using a complex scalar field 
theory, superfluidity arising via a Bose-Einstein condensation. Note, however, that a 
fundamental description would treat 4 He as a tightly-bound state of fermionic con- 
stituents, not too dissimilar in spirit to Two-Color QCD. 3 He by contrast is a fermion, 
and superfluidity in this case is believed to arise via a BCS instability resulting in 
a condensation of weakly-bound Cooper pairs. We might thus consider superfluidity 
in the NJL model as a relativistic generalisation of this phenomenon. It is impor- 
tant to note, though, that due to short distance repulsion between helium atoms the 
BCS wavefunction in 3 He is actually p-wave, resulting in ground states described by 
a complicated order parameter and many interesting topological excitations ||14|| . In 
the NJL model studied here, the corresponding qq wavefunction is a scalar s-wave, 
and any superfluid might be expected to behave more like 4 He. 

Superfluidity in a relativistic model similar to the NJL model has been studied 
using mean field techniques in || . In our previous work [15|, [U| we have attempted to 
identify superfluidity in the 2+ld NJL model using non-perturbative numerical lat- 
tice simulations. Apart from the obvious computational gain, we chose this particular 
dimensionality because the model has a non-trivial continuum limit |17| , There- 
fore, in contrast to effective descriptions such as the Landau-Ginzburg theory, the 
condensed matter described in this approach is formed from the elementary quanta 
of an interacting field theory. Our results have not supported the expected scenario 
outlined above, raising the question of whether important physics is neglected in the 
self-consistent approach. Although there is evidence for enhanced diquark pairing in 
the scalar isosinglet channel in |15| , we have not succeeded in finding an unambiguous 
signal for a condensate (qq) ^ 0. Rather, in the high density phase the condensate 
appears to vanish non-analytically as a function of diquark source strength j, sug- 
gesting critical behaviour [Q. Studies of the excitation spectrum in the spin- 1 ; sector 
reveal a sharp Fermi surface and no evidence for a gap A ^ [|T^j . The purpose of the 
present paper is to present these results in greater depth, and to attempt to interpret 
them. As we shall see, it may be possible to attribute the unconventional signals to 
the specifically two-dimensional nature of the system being studied, which thus bears 



many resemblances to superfluidity observed in thin helium films [|19| . We will argue 
that neither long range order (qq) ^ nor a gap A / are necessary attributes of 
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a superfluid. Instead, the critical behaviour observed in [16] results from long range 



coherence in the phase of the diquark wavefunction, which appears to be a sufficient 
condition for thin film superfluidity [POR . 

In Sec. [| we review the formulation and numerical simulation of the lattice NJL 
model in 2+ld with non-zero chemical potential /i, paying particular attention to 
the introduction of diquark source terms jqq via the use of a Gor'kov basis 0. It is 
thus possible to define diquark observables which are measurable on a finite system; 
Sec. |3] reviews numerical results for the diquark condensate (qq{j)), the associated 
susceptibilites, and diquark masses. Critical behaviour in the high density phase 
tiB > is identified after extrapolating results for the first two quantities to the 
zero-temperature (ie. L t — > oo) limit, leading to consistent estimates for the critical 
exponent conventionally denoted S which vary with \x. In Sec. f| this behaviour is 
discussed in analogy with that of the 2d XY model, in which long range order is 
washed out by spin wave excitations and which also displays critical behaviour in a 
continuous parameter range. It is argued that in such circumstances persistent flow, 
the defining property of a superfluid, can only be disrupted by excitations costing 
infinite energy in the thermodynamic limit. Sec. |5| presents the results of a study of 
the dispersion relation E{k) of spin-| quasiparticle excitations in the dense phase, 
revealing a sharp Fermi surface for the first time using lattice methods. There is 
no evidence for particle-hole mixing or a non-vanishing gap as the source strength 
j — > 0. Instead, the results in this sector are consistent with a normal Fermi liquid 



of the type first discussed by Landau [ElL 22]; in particular it is possible to estimate 



both Fermi momentum kp and velocity (3p as functions of /x and to show that these 
depart from their free-field values, yielding information on quasiparticle interactions. 
Conclusions and suggestions for further work are outlined in Sec. |6] 

2 The Lattice Model 

2.1 Formulation and Symmetries 

The model studied in this paper is a lattice transcription of the NJL model in 2+1 



dimensions, identical to that studied in [15, 16]. It is defined by the Euclidean action 



S — Sf er + St, os : 

S fer = J2xMmX + 3X tr T2X + 3XT 2 x tr ; S bos = 4E tr$t$ > (2- 1 ) 

X 9 x 
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where x> X are isospinor fermionic fields defined on the sites x of a 2+ld lattice, and 
$ = cr + i-K.T is a 2 x 2 matrix of bosonic auxiliary fields living on the dual sites x. The 
kinetic operator M has the standard form for staggered lattice fermions interacting 



with scalar fields 118 



Mm] = \t pq 



2 

+ ls xy J2[a(x)5^ + ie(x)7r(x).r^). (2.2) 

(x,x) 



( eM <Wo _e ^ 5 yx-o)+ Vu(x)(5 yx+0 -5 yx ^) + 2m6 xy 

u=l,2 



The parameters are bare fermion mass m, baryon chemical potential fi, and coupling 
g 2 . The r are Pauli matrices normalised to tr^Tj) = 2<5^ acting on isopin indices 
p,q = 1,2. The symbols rj v (x) denote the phases (—l} x o+-+x v -i ; e ( x ) the phase 
(— i^ x o+xi+x 2 ^ anc [ {x,x) the set of 8 dual sites neighbouring x. Integration over the 
auxiliary $ fields leads to an equivalent action in terms of fermions which self-interact 
via a four-point contact term proportional to g 2 , corresponding to the interaction of 



the NJL model 23 



In addition to the usual NJL interactions, Eqn. (2T) contains diquark and anti- 
diquark terms proportional to source strengths j and J respectively. These have been 
introduced to enable the measurement of the diquark condensate (x trr 2X) on a finite 
system, in precise analogy to the role of the bare mass m in the measurement of the 
chiral condensate (xx)- To proceed, we define the bispinor ^> tr = (x tr , x)> an d rewrite 
the fermion action as a quadratic form Sf er = ty tr A^, where in this Gor'kov basis 
the antisymmetric matrix A is 

The fermion fields may then be integrated out to yield the following Euclidean path 
integral: 

Z = J DaDn Pf(2^[$, j, j]) exp -S 6os [$], (2.4) 

where the pfaffian Pf(Q) = \/detQ. Note that this expression differs from the (in- 
correct) version given in by a physically irrelevant factor of two; for convenience 
we will stick with the current notation, but note that if the source is interpreted as a 
Majorana mass A, then j = A/2 ||. 

The model described by the action fl2.iy2.2j ) has an SU(2) L <g> SU(2) R ® U(1) B 
global symmetry. Defining projection operators V e / = \{\ ±e(x)) onto even and odd 
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sublattices respectively, we have 



X 



{V e U + V V)x; X ^ x(VeV^ + V U^); $ i— > V<&U^ [U,V E SU(2)] : (2.5) 



X^e^x; X^Xe~ ia [e ta E U(1) B ]. (2.6) 

The symmetry ( |2.5|) is broken to a diagonal SU(2)y of isospin with U = V in (|2.5|) , 
either explicitly by a bare fermion mass m ^ 0, or spontaneously by the generation 
of a chiral condensate (xx) 7^ by the model's dynamics. For fi = this occurs for 
a sufficiently strong coupling g 2 > g c ~ 1.0a, where a is the physical lattice spacing 
||18||. In the chirally broken phase the fermions have a dynamically generated mass 



E ~ (a) = 2"(xx). Since Ea — > as g 2 — > g 2 , a continuum limit may be taken at this 
critical point. A remarkable feature of the 2+ Id NJL model is that the continuum 
theory so obtained remains interacting |]I7|, |TE| . As in our previous studies [|T5|, [T(J , the 



simulations in this paper were performed with g 2 = 2.0 corresponding to Ea = 0.71, 
implying that we are rather far from the continuum limit. 

For \x ^ the model is known to exhibit a strong first-order transition to a chirally 



symmetric phase [23|; for our realisation this occurs at a critical fi c ~ E ~ 0.65, as 
shown in Fig [I] ||15|| . Chiral symmetry restoration is accompanied at this point by the 
onset of a non-vanishing density of baryon charge in the ground state, signalled by a 
condensate 

n B = = lixix^xix + 6) + x(x)e^x(x - 6)) > 0. (2.7) 

Existing numerical evidence suggests these two a priori distinct transitions are coinci- 
dent f25|. For /i > /i c , E ~ 0, and the density follows the approximate form n B oc /i 2 , 



the behaviour expected for massless states populating a 2 dimensional Fermi sphere 
of radius E F = fi, until it gets close to its saturation value of one quark of each isospin 
per lattice site. For our realisation, n# ~ 0.25 quarks of each isospin label per site at 
/i = 0.8, corresponding to a physical density of about 80 fm -2 assuming a constituent 
quark mass E of 300 MeV; by /i = 0.9 this has risen to % ~ 0.4, corresponding to 
125 fm" 2 . 

The question which will occupy us in this paper is the nature of the high density 
phase present for ji > /i c , and in particular whether the U(l)s symmetry (|2.6j ) is 
spontaneously broken by the generation of a diquark condensate which we will gener- 
ically denote by (qq) ^ 0. In a BCS condensation, the participating diquark pairs 
come from the neighbourhood of the Fermi surface. The resulting ground state is 
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Figure 1: Chiral condensate £ and baryon charge density % asa function of fj, for 
a 16 2 x 24 system with g 2 = 2.0, m = 0.01, j = 0. 

separated from excited states by an energy gap A analogous to the constituent quark 
mass S in a chirally-broken vacuum; mean-field calculations of a phenomenologically 
inspired 3+ld NJL-type model predict that for quark matter A is of the same order 
as E 0. As well as being massive, the quasiparticle excitations above the ground 
state carry indefinite baryon charge due to the U(l)s breaking. Physically, this means 
that the quasiparticle is a coherent superposition of particle and hole states. In Ref. 
P~5| ] diquark timeslice correlators Y^x(QQ(^j Q)qq(%i 0) m var i° us plausible condensa- 
tion channels were studied for \l > fi c , and evidence for pairing found, in the form of a 
plateau whose height did not decrease with Euclidean time separation t, for the scalar 
SU(2)i®SU(2)^j singlet channel qq = x trr 2X- However, the naive interpretation of 
diquark condensation via the clustering hypothesis, namely that 

}im(qq(0)qq(t)) = \(qq)\ 2 , (2.8) 

t — >oo 

was ruled out because the plateau height did not scale in the expected way, ie. exten- 
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sively in the spatial volume L 2 S . To clarify the situation, Ref. |TBj introduced diquark 
source terms, as in fl2.1|) , making direct measurements of (qq) possible; we now review 
the 'standard' signals which might be expected if diquark condensation occurs. 
Firstly we define diquark operators qq± via 



qq±(x) 



X r (x)T 2 x{x) ± x(x)t 2 x r (x) 



(2.9) 



with corresponding source strengths j± — j ± J. It is readily checked that qq± are 
invariant under SU(2)i®SU(2)# rotations ( ^.5[ ), but rotate into each other under 
U(l)s (2J3). In terms of 4-component spinors ip, the operators (|2.9| ) may be written 

m 



V r (C<y 5 ® r 2 ® r 2 )^ ± ^75 ® r 2 ® r 2 )^ 



tv 



(2.10) 



where the first matrix in the tensor product acts on spinor indices, the second on a 
2-component flavor structure which is implicit in the staggered fermion approach [^6 



and the third on the explicit isospin index introduced in (|2.1| ). The charge conjugation 
matrix C is defined by C r ~f fJ- C^ 1 = — 7*. The diquark condensate is now given by 

1 d\nZ 



-) 



AV 



(2.11) 



v du 

and is calculable using standard lattice techniques, such as the use of a stochastic 
estimator for the diagonal elements of the inverse matrix. The non-vanishing of 
( [2. 11| ) in the limit j + — > is a criterion for the spontaneous breakdown of U(l)# 
symmetry. Furthermore, if we define susceptibilities 



X± = X^±(%tf±( x ))> 



(2.12) 



then it is straightforward to derive a Ward identity analogous to the axial Ward 
identity for the pion propagator: 

(??+) 



X- 



J+ 



(2.13) 



On the assumption that the dominant contribution to X- is from a simple pole, then 
gg_ couples to a Goldstone mode whose mass M_ vanishes in the zero source limit as 
If we similarly attribute x+ t° a "Higgs mode", then the ratio x+/X- provides 
an alternative means of distinguishing possible symmetry-breaking scenarios in the 
limit j + —>■ 0: 

X+ _ ( 1, if U(1) B manifest; 
^7 10, if U(1) B broken. 



R 



lim • 



(2.14) 
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In Sec. [3] we will present numerical results for these quantities and discuss to what 
extent the above considerations help in describing the high density phase of the 2+ld 
NJL model. 



2.2 The Simulation 



Numerical simulation of the path integral ( |2.4j) requires some discussion of how to 
deal with the pfaffian. First let us find the condition that det2*4 is real. Using the 
property of a block square matrix 



det 



fX Y 

\w z 



\ = detXdet(Z - WX~ l Y), (2.15) 

and the property r 2 Mr 2 = M* which follows from fl2.2p , we deduce 

det2^l = det(4jJ+M t M), (2.16) 

and is hence real and positive if jj is chosen real and positive. In all our work we 
choose j = j real which satisfies this condition. It follows that Pf(2^4) is real. In 
fact, one can go further and argue that it is also positive as follows. In the limit 
j,j — > 0, Pf(2^4) reduces to detM, which can be proven both real and positive using 
an argument, identical to that used for SU(2) lattice gauge theory with staggered 
quarks in the fundamental representation ||, showing that any complex eigenvalue of 
M is accompanied in the spectrum by its conjugate, and any purely real eigenvalue 
is doubly degenerate. Relation ( |2.16|) , however, shows that det2^4. can only increase 



and hence cannot change sign once jj > 0; it follows that Pf(2^4) can be consistently 
chosen real and positive ||. 

Despite this reassuring property, in our simulation we chose to use oc 
Pf 2 (2^l) as the measure, corresponding to two staggered lattice fermion species, for 
consistency with the model of |l5j which is recovered in the limit j — ► 0. It can 
be shown that in the continuum limit the model contains Nf = 4 species of four- 



component Dirac fermions |26j . The simulation is performed using a hybrid molecular 
dynamics "R" algorithm [27]], in which the square root is taken by inserting a factor 
of \ in the force term derived from a local action; note that we were able to debug 
and tune the code by also checking against an exact algorithm for the case Nf = 8. 
In all cases we used a molecular dynamics timestep At = 0.04, and never saw any 
evidence of departure from equipartition of energy. We performed simulations on 
lattice sizes L 2 S x L t = 16 3 , 16 2 x 24, 24 3 , 32 3 , and in one case 48 3 , with the coupling 
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g 2 fixed to 2.0 as described above, and bare Dirac mass m fixed to 0.01 to assist with 
the identification of chirally broken and restored phases. A typical run is over 0(400) 
HMD time units with mean refreshment interval 1.0; data were taken every two units. 
The cost of the simulation rises considerably in the chirally-restored phase where the 
diagonal elements of M, proportional to (a), are small, particularly as j — > 0. The 
48 3 point at /i = 0.8, j = 0.025 required approximately 40 SGI Origin2000 processor 
days. 

In Sec. p.l| we review the behaviour of the model as a function of fi, and present 
results for (qq) taken in the "partially quenched" approximation in which j ^ only in 
the measurement, and the simulation performed using an exact algorithm with j = 0. 



Our studies with full pfaffian dynamics, presented in Sees. p72] , p73] and p^4| , focussed on 
two representative points in the chirally symmetric low density phase at /i = 0.0, 0.2, 
and on two values in the high density chirally symmetric phase fi = 0.8, 0.9. We used 
j = J ranging in value from 0.3 down to 0.025. In our studies of the quasiparticle 
spectrum discussed in Sec. |5| we performed runs on 32 3 at 4 additional values of 
fi e (0.8,0.9). 



3 Numerical Results for Diquark Observables 
3.1 Partially Quenched Results 

In this section we discuss the direct numerical measurement of diquark condensates 
(qq), which the Gor'kov basis discussed in the previous section makes possible. To 
warm up we consider the partially quenched approximation, in which the source 
strength j is set to zero in the updating of the {$} configuration, which thus pro- 
ceeds via an exact Hybrid Monte Carlo algorithm as in |I5| , but is non-zero in the 
measurement routine so that (qq) can be measured via ( |2.11|) . Since at a given /i a 
single simulation serves for all j, this approach is fairly cheap and hence good cover- 
age of the /i axis is practicable. Our results are shown in Fig. 0, where open symbols 
denote data taken in the low density phase \i < /i c , and filled symbols are from the 
dense phase. For n < /i c (qq+) varies approximately linearly with j, implying that a 
smooth extrapolation to the origin is possible, and hence the condensate vanishes in 
the j — > limit. A striking jump occurs between \i = 0.65 and ji = 0.7, and for values 
of fi in the dense phase (qq(j)) is markedly more curved. This behaviour cannot be 
taken as evidence of diquark condensation, however; one should expect discontinuities 
in all physical observables on different sides of a first order phase transition. Despite 
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Figure 2: Partially quenched results for the diquark condensate (x* r r 2 x) as a function 
of source strength j for various values of /i; unless otherwise shown, data was taken 
on a 16 2 x 24 system. 

the curvature of the lines a smooth extrapolation to the origin consistent with unbro- 
ken baryon number symmetry at high density is still plausible. Another possibility, 
suggested by the \x = 0.8 data from a 30 3 lattice shown in Fig. is that the sym- 
metry breaking is masked by a finite- volume suppression as j ' ^ 0. To explore the 
behaviour for fi > \i c in more depth data from lattices with several distinct L s and 
L t are needed. We chose to perform this using "unitary" data generated using the 
correct measure (|2.4j ) for a limited number of different fi, as described next. 

3.2 Diquark Condensate 

In Fig. |3] we show (qq+) data taken from simulations using the full pfaflian measure 
( |2.4|) at two values of fi from each phase on a 32 3 lattice. The results resemble those of 
the partially quenched approach shown in Fig. ^| both qualitatively and quantitatively. 
The curvature of (qq+(j)) in the dense phase seems to become more pronounced with 
increased fi, to the extent that by j ~ 0.2 the results at \i = 0.9 actually undershoot 
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Figure 3: Results for (qq+) vs. j for various values of /i from a full simulation on a 
32 3 system. 
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Figure 4: Extrapolation of (qq+) at /i = 0.9 to the thermodynamic limit. 
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those at fj, — 0.8. As remarked in the previous section, there are significant finite 
volume effects in this phase. Fig. |] shows fi = 0.9 data from simulations on 16 3 , 
24 3 and 32 3 lattices. The equivalent data for fj, = 0.8, including a single point from 
48 3 , is tabulated in Table [l] and plotted as Fig. 2 of [|T6 1 . Empirically, we find by 



comparing data from 16 3 , 16 2 x 24 and 24 3 lattices that the dominant correction on 
a L 2 x L t system appears due to finite L t , suggesting a specifically thermal origin. 
This motivates an extrapolation to the thermodynamic limit which is linear in 1/L t ; 
at smaller j, however, the data depart significantly from this trend. 

Assuming a 1/Lt scaling, we extrapolated the data from 24 3 and 32 3 lattices to 
estimate (qq + (j)) in the thermodynamic limit. The results are shown on a log-log 
plot in Fig. [|, together with unextrapolated data from the chirally broken phase at 
\x = 0.0, 0.2. Remarkably, there is a reasonably wide interval j G [0.05, 0.2] within 
which the plot is approximately linear, indicative of a power-law scaling 

(qq + (j))<xj a . (3-1) 

Fits to (£]]) in this range yield a = 0.314(2) for // = 0.8, with x 2 /dof = 2.3, and 
a = 0.213(3) for \i = 0.9, with x 2 /dof = 0.4. In both cases this is clearly distinct 
from the linear (ie. a = 1) behaviour observed at low density. For j outside the 
fitted range, the data start to fall below the fitted line; we ascribe this to scaling 
violations for j > 0.25, as perhaps revealed by the crossing of the curves in Fig. [| 
and for j ^ 0.035 to non-thermal finite volume effects, eg. due to an insufficiently 
large explicit Majorana mass, as perhaps indicated by the different systematics of the 
16 3 point in Fig. ^. Unfortunately our resources have not permitted further systematic 
study of this point. 
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Figure 5: hx(qq + ) vs. Inj, showing evidence for power-law scaling in the dense phase. 



Assuming the validity of the form ( |3.1|) , we draw two conclusions. Firstly, the 
non-analytic behaviour is reminiscent of the power-law scaling observed at a critical 
point of a thermodynamic system. For a spin system at its critical temperature the 



spontaneous magnetisation Ai scales with applied magnetic field h as M. oc hs [^8 
for a fermionic model exhibiting chiral symmetry breaking the equivalent relation is 
{ipip) oc mh [jl^l . We thus identify critical scaling, with 5 = a' 1 . Secondly, Fig. |2] 



leads us to expect that critical behaviour is generic in the dense phase, but with the 
exponent 5 varying continuously with chemical potential /i, taking the value 5 ~ 3 at 
fi = 0.8 and ~ 5 at /i = 0.9. This suggests a line of critical points for /i > /i c . The 
origins of such behaviour and its consequences for superfluidity will be elaborated in 
Sec. [|. In an attempt to find further evidence for criticality, however, we now switch 
attention from one- to two-point functions in a study of the various susceptibilities. 
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3.3 Susceptibilities 

Next we examine the diquark susceptibilities x± denned by ( |2.12|) , which may be 
expanded to 

x 

Mx tr T2X(0)XT2X tr (x) + XT2X tr (0)x tr T 2 x(x)}. (3.2) 

A generic susceptibility may be expressed as the sum of two connected contributions 
corresponding to the two possible Wick contractions, 



X 



<(trr<^) 2 > - (trreu 2 + (trSo^r) = x s + x ns , (3.3) 



where Q = A^ 1 is the Gor'kov propagator and T projects out the appropriate compo- 
nents. By analogy with meson physics we label these contributions "singlet" and 
"non-singlet" respectively. Estimates for x± are made with the same stochastic 
method used for (qq+), and are plotted for /i = 0.8 on a 32 3 lattice in Fig. [| Apart 
from the observation that X- ~ X+i no other trend is apparent in the data, which 
are noisy and quite possibly consistent with zero. In the following we ignore x± an d 
assume x± — X±- This is in marked contrast with the behaviour observed in studies 
of chiral symmetry breaking in 2+ld fermionic models, where singlet contributions 
to the relevant susceptibility are significant [29| or even dominant [|10|j . 



Restricting attention to the non-singlet pieces, it is not hard to show using the 
properties of Q reviewed in Sec. |5| below, that the first expectation value on the right 



hand side of (|3.2|) is negative and vanishes in the limit j — > 0, whereas the second is 
positive and in fact corresponds to the diquark correlator examined for j = in [|15|| . 
We conclude that \x~\ > \x+\- Data for x± from a 32 3 lattice are shown in Fig. [7[ 
Data at /i = 0.0 are indistinguishable from those at /x = 0.2 on the scale plotted. 
Note the difference of scale on the vertical axis between Figs. ^ and 0. We should also 
comment that the x~ s data when checked against (qq + ) saturate the Ward identity 



( [2.13] ) within errors. Both observations justify our neglect of x±- 

In a conventional symmetry breaking scenario X- should diverge in the thermo- 
dynamic and j — > limits according to (|2.13j ) , whereas x+ could in principle remain 



finite. The ratio R defined in ( |2.14| ) is expected to vanish as j — > if U(l)# symmetry 



is spontaneously broken by the ground state, and to approach unity if the symmetry 
remains manifest. In order to investigate this it is once again necessary to take ac- 
count of finite volume effects. There is no appreciable effect for ji < fi c , but in the 
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Figure 8: Susceptibility ratio R at fi = 0.8 for various lattice sizes. 

dense phase the variation with lattice size is considerable, as shown in Fig. [8]. Once 
again, an extrapolation oc L^ 1 seems plausible, and indeed in this case a linear fit 
to the data from all three available volumes proved acceptable. The accumulation of 
the resulting intercepts in the region R ~ 0.3 is striking. 

The extrapolated results for R as a function of j are shown in Fig. H In the chirally 
broken phase the results support R tending smoothly to one as j — > 0, consistent with 
unbroken baryon number symmetry. The behaviour in the high-density phase is very 
different; the accumulation of intercepts in Fig. |8| manifests itself as a plateau for 
j ^ 0.075. For smaller j the ratio shoots sharply upward towards one. This can 
be attributed to a finite volume artifact, since we know that terms in ( |3.2| ) of the 
form (qq{0)qq(x)) which split the degeneracy between x+ an d X- necessarily vanish 
as j — > away from the thermodynamic limit. In this regard it is encouraging that 
these non-thermal effects manifest themselves in the same range j ^ 0.05 observed 
for the condensate measurements of the previous section. We are thus motivated to 
attempt a linear extrapolation to j = for the data with j G [0.75, 0.2]. The fits are 
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Figure 9: Susceptibility ratios R extrapolated to infinite volume for variuos \i = 0.8. 

of excellent quality and yield R{j = 0) = 0.29(2) for \l = 0.8 and R(j = 0) = 0.17(1) 
for fi = 0.9. 

Measurements of the diquark condensate of Sec. |3]^ support a power-law form 
(qq+) oc j a (p.l|). If this is the case, then the relation x+ = d(qq + ) /dj + together with 
the Ward identity ( ^l3|) imply |3 



<91n gg + 

^0) = "^j — : — = (3-4) 

which crucially is independent of j. The plateaux of Fig. ||] clearly support this 
interpretation; moreover the values we obtain for R(j = 0) are in surprisingly good 
agreement with those from fits to ( |3.1| ). The susceptibility measurements thus provide 
an independent corroboration of the hypothesis that the system is critical for > fi c . 



19 



3.4 Diquark Masses 

Our final numerical study in this sector is of the spatial behaviour of the diquark 
correlation functions, in an attempt to estimate the masses M± of diquark bound 
states. For brevity we only consider \i > fi c , in which case M± are probably best 
thought of as the energies required to excite a diquark pair above the ground state. 
We have restricted our attention to the zero-momentum timeslice correlator C±(t) = 
Sx(99±(0, 0)qq±(x, t)), so that the excited state must consist of quarks with equal and 
opposite momentum k. Recall that in the presence of a Fermi surface, only quarks 
with k ~ hp can be excited; the measurements presented here are not sensitive to 
this restriction, although it will prove a decisive factor in the quasiparticle study of 
Sec. [5|. As in the previous section, we ignore "singlet" diagrams in calculating C±. 



Fig. [10] shows the correlators for j = 0.025 and 0.25 at fx = 0.8. By virtue of its 



definition, C± is clearly symmetric under time-reversal, in contrast to the correlators 



studied in |15|]. It is also clear, as expected, that |C_| > \C + \, and that the difference 
between them grows with j. Close inspection of Fig. |10] suggests that a standard 
simple-pole fit to C±(t) will not succeed unless a constant term is included; we have 
therefore attempted fits of the form [pjj 

C±(t) = P ± (exp(-M ± t)+exp(-M±(L t -t))^+Q ± . (3.5) 

The plateau height Q + would by the cluster property be proportional to | (qq+) | 2 if the 
condensate formed; however, the analysis of |L5] showed that at j = 0, Q + does not 
display the required extensive scaling with two-dimensional spatial volume. There is 
no obvious theoretical interpretation for Q_. 

Fig. [XT] shows M±(j) for li = 0.8,0.9 resulting from fits of ( |3.5| ) to timeslices 



5-26. In most cases the x 2 /dof was ^2.0 and in no case exceeded 6.0. M ± is 
found to increase almost linearly with j, maintaining a roughly constant difference 
M + — M_ ~ 0.08 for j > 0.05. For smaller j the curvature in the plots suggests the 
two states may become degenerate as j — > 0. A linear extrapolation to j = yields 
M + (lm = 0.8) ~ 0.23, M + (/i = 0.9) ~ 0.21, values of the same order of magnitude but 
slightly lower than those obtained directly at j = on a 16 2 x 24 lattice in JOj) (note 
that the symbols M± have a different meaning in that paper). The fitted values for 
Q± vary considerably over the range of j explored: eg. for li = 0.8, Q + rises from 
0.102(1) x 10" 5 at j = 0.25 to 0.315(5) x 10" 1 at j = 0.025; in the same range Q_ 
rises even more dramatically from 0.49(6) x 10~ 7 to 0.279(7) x 10 _1 . 
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Figure 11: Diquark masses extracted from fits to ( |3.5| ) on a 32 3 lattice. 
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The most important feature of Fig. [11] is that there is no evidence for M_ vanishing 
as j —> 0, as might be expected if gg_ coupled to a Goldstone mode as a result of 
broken U(1)b symmetry. One might argue that the ad hoc inclusion of Q_ in the fit 
( |3.5|) results in artificially high values of M_ ; in any case, the conclusion remains that 
simple pole fits to C-(t) corresponding to a weakly interacting Goldstone boson in 
this channel fail drastically. The scaling form ( |3.1| ) combined with the Ward identity 
( f2.13j ) imples a massless degree of freedom as j — > 0, and hence long range correlations, 
in both qq_ and (since symmetry is restored in this limit) qq + channels; they must, 
however, be strongly interacting and hence short-lived states. 



4 Criticality and Superfluidity 

Having established that in the limit of vanishing source there is no diquark conden- 
sation at high density, but instead a critical phase with the scaling of the condensate 
with the source governed by an exponent 6 varying continuously with /z, we now dis- 
cuss the implications for possible superfluid behaviour of the 2+ld NJL model. In 
fact, this result is in accordance with well-known theorems that long range ordering 



of a two-dimensional system with a continuous global symmetry is impossible [|3l 



In the current context a particularly appropriate statement of the theorem is due 



to Hohenberg [p2 |, who explicitly considers the case of a composite order parameter 
via Cooper pairing in a low- dimensional fermion superfluid. Long-wavelength fluctu- 
ations of the phase 9 of the would-be condensate always wash out the order in the 
zero source limit. In field theoretic language, in two dimensions infra-red divergences 
dictate that the Goldstone pole in the transverse susceptibility predicted by a naive 



application of fl2.13| ) is replaced by a softer divergence consistent with a power-law 



decay of the correlator 

\im(qq(0)qq(r)} oc (e !9(0) e" tf(r) ) oc — , (4.1) 

where r\ is another critical exponent, implying a massless but strongly- interacting 
mode and long-ranged phase correlations. Note that direct numerical tests of (|4.1| ) 
would require data from spatial diquark correlators, in contrast to the temporal cor- 
relators explored in Sec. |3.4| . There are also technical difficulties in taking the limit 

The best known example of a system with a critical phase is the 2d 0(2) spin 
or XY model, which is similar in that long range order would also spontaneously 



22 



break a U(l) global symmetry. The critical behaviour occurs for T < Tbkt-, the 



temperature of the celebrated Berezinskii-Kosterlitz-Thouless transition [£0|, [33[] . The 
physical picture can be explained as follows; on the assumption that the interaction 
strength is a periodic function of the difference in angle 9 between adjacent spins, 
and is approximately gaussian in neighbourhood of its minima, then an effective 
Hamiltonian can be written as 



H XY [6,m\ = jJ2(d^0(x)) 2 -27rJ^m(x)ln 

Xfl x,y 



x-y 



r 



m(y), (4.2) 



J being the nearest neighbour coupling. In addition to the 8s the Hamiltonian depends 
on integer-charged vortices m(x) located on the sites of the dual lattice. The vortices 
are topological excitations of the spins which interact via a Coulomb potential which 
is logarithmic in two dimensions, ensuring that all configurations with finite Hxy are 
overall charge-neutral, ie. contain as many anti-vortices as vortices. The parameter 
r is the "core size" of the vortex, which can be considered of the same order as the 
lattice spacing. Now, at low temperatures the second term of ( (4.2| ) strongly suppresses 
well-separated vortex - anti-vortex pairs, and the model's dynamics are dominated 
by small-amplitude fluctuations of 8, the so-called spin waves. Phase correlations are 
governed by ( |4.1| ) with t](T) = T/4ttJ, implying a critical phase with continuously 
varying exponent. At the critical temperature Tbkt, the vortex entropy begins to 
dominate the free energy of the system, and vortex pairs of arbitrary separation 
form. The resulting vortex plasma screens long range correlations resulting in a finite 
correlation length for T > Tbkt- 

Next we discuss the relation with superfluidity. We can rewrite the diquark op- 
erator qq + (x) = 4>{x) = 4>oe %e ( x \ where the constant 0o is the density of quark pairs 
participating in the condensate and 9 the local phase of the diquark operator. In this 
form qq + (x) can be regarded as a bosonic macroscopic wavefunction for the condensed 
pairs. We now identify a superfluid current J Sfl via 

J a(t oc -~ (0*<9 M - (<9 M 0*)0) = K s d,6. (4.3) 

The constant K s must be determined empirically. In the non-relativistic limit it is 
given by 

K s = ^n s (4.4) 

where M is the mass of the current- carrying atomic species (M( 4 He) or 2M( 3 He) 
in the case of the two known superfluids), and n s a parameter called the superfluid 
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density, which for an interacting system need not coincide with the charge density of 



the particles in the condensate p2| . In turn this enables the definition of a superfluid 
velocity v s = ^W. For a relativistic system the relation v s = ^V6* can be shown 
to hold for diquark pairs for small v s f34fl . 

Now in the static limit, relation ( f4.3|) implies that the flow is irrotational, viz. 
V x J s = 0, and hence the circulation k = § J s .dl around any closed path vanishes 
unless either the condensate is somewhere singular within the contour, ie. 0o — 0, or 
the space is non-simply connected. In either case the requirement that 9 be single- 
valued results in the quantisation of circulation: n = 27rnK s , with n integer. In the 
case of a singularity in <p the physical realisation of ft 7^ is a vortex, with a non-zero 
radius tq within which the normal phase is restored. Superfluid vortices experience 
long-ranged mutual interactions; for a two-dimensional system such vortices can be 
identified with the vortices of the XY model, and are expected to be governed by 
an effective Hamiltonian of the same form as ( [4.21) . An example of a non-simply 
connected space would be a finite system of dimension L x x L y with periodic boundary 
conditions; in this case k^O implies a uniform supercurrent 



J s (n x ,n y ) = 2ttK s + ^yj . (4.5) 

The crucial point is that the resulting flow patterns are topologically stable, implying 
the system's energy must be greatly increased to change k ||20|| . For instance, in 
order to change n x by one unit, a vortex - anti- vortex pair must be created and 
the vortex moved in the ^-direction right around the system before being allowed to 
reannihilate with the anti-vortex. In so doing the system must be brought through a 
saddle-point configuration in which the pair is separated by half the system extent; 
from ( |4.2| ) the energy required ~ 27rJln(Lj / /2r ). Since this diverges with the size of 
the system, we infer that the circulation is stable and hence the current J s persistent, 
implying superfluidity. We conclude that the critical phase of the XY model, and 
by extension critical behaviour in any two dimensional system with a U(l) global 
symmetry, describes superfluidity despite the absence of a condensate. Long range 
order of the phase 9 is not necessary; phase coherence, as expressed by (|4.1| ), is 
sufficient. It is noteworthy in this regard that some of the most precise tests of 
the universal predictions of the XY model have come from studies of thin films of 
superfluid 4 He |19| . 



Since we have used universal features of vortices and spin waves to argue that 
critical behaviour implies superfluidity in two dimensions, to justify the application 
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of these ideas to the NJL model we should address the issue of why our 5 is not 
consistent with that of the XY model as revealed by a renormalisation group analysis 



35] , which predicts 

<5>15 ; 77 <i (4.6) 

with equality holding as T — > Tbkt-- First we note that dimensional reduction, which 
predicts that the critical thermal properties of 2+ld systems should be governed by 
the 2d spin model with equivalent global symmetry, does not apply in this case (see 
PEfl for a recent numerical study of a four-fermi model with U(l) axial symmetry at 
T / which does appear consistent with the BKT scenario). Rather, the feature 
which permits us to use a two-dimensional effective model is the static nature of the 
phase fluctuations, ie. d t 6 ~ 0, as evidenced by the plateaux observed in the large-t 



behaviour of diquark correlation functions in [15| and Sec. |3.4| . Note in addition that 
we have needed the limit L t — > 00 rather than L t — > 0. Since the number of accessible 
Matsubara modes in our simulations remains large, the fermions need not decouple 
(indeed, there remain light fermionic excitations, as we shall see in Sec. [5]) and we 
should not expect the model's dynamics to be described by a purely bosonic effective 
action. Symmetry breaking via a composite order parameter is qualitatively different 
from cases where the order parameter is an elementary field; there is a wealth of 
evidence, both analytical and numerical, that bosonic and fermionic models with the 
same patterns of global symmetry breaking belong to separate universality classes in 
dimensions up to [17], [18], |29], [37J and even including |23|] four. 

Our simulations have yielded information only about the critical exponent 5 in 
the dense NJL model. Our results 5(p, = 0.8) ~ 3, S(fi = 0.9) ~ 5 are consistent 
with a critical phase for fi > fi c . The lower numerical values as compared to (4.6) 
typify the distinct nature of symmetry breaking via a composite order parameter as 
discussed above. Note, however, that although S decreases with /x, the analagous 
I^bkt at which superfluid vortices unbind may not be physically accessible; most 
probably chiral symmetry breaking at /i = [l c happens first. Finally, although we 
have not yet found a method to measure 77, it is interesting to estimate its value using 
the hyperscaling relation 5 = (d + 2 — 77) / (d — 2 + rj) [|28] . Since we have assumed an 



effective dimension d = 2 for the critical dynamics, the appropriate relation is 

S = (4.7) 
V 

yielding T](fi = 0.8) ~ 1 and T](fi = 0.9) ~ 0.7. Note that had we used an effective 
dimension d = 3, the prediction for 77 at \i = 0.9 would be almost vanishing. 
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5 The Quasiparticle Spectrum 



In this section we study the spin-| sector by examining the Gor'kov propagator 
Q = A^ 1 as a function of fi and j. For fi < \x c the fermion excitations are sim- 



ply related to those at // = 0, as reported in [fLl}| . For /i > /i c , however, the ground 
state of the model changes radically, and is characterised by a Fermi surface with 
energy Ep and characteristic momentum hp. A generic description is Fermi liquid 
theory pl| , p2| , in which excitations with momentum such that |A; — &p| <S are 
quasiparticles whose mass need not coincide with that of the fundamental quanta. If 
a BCS condensation occurs, then the lowest energy excitation may be separated from 
zero by a gap A, and the quasiparticles may not be eigenstates of baryon number, but 
instead some kind of particle-hole superposition. Analysis of Q using standard lattice 
spectroscopy techniques yields information on the quasiparticle dispersion relation 
E(k), thus probing A and, more generally, the nature of the model's Fermi surface. 

We begin by making some general observations about the Gor'kov propagator. 
Write 



\N(x,y) A(x,y) 

where each element denotes a 2 x 2 matrix in isospace. Our notation signifies that the 
propagator contains both "normal" (q(x)q(y)) and "anomalous" (q(x)q(y)) compo- 
nents, together with their barred counterparts. On a finite system, A and A vanish in 
the limit j — > 0, and N, N become proportional to the usual fermion and anti- fermion 
propagators. The number of independent components of Q is constrained by certain 
identities. For instance, N is proportional to an element of SU(2), implying N22 = N* x 
and N21 = —N* 2 , with similar relations for N. In the anomalous sector, however, it 
is t 2 A which resembles an SU(2) matrix so the equivalent identities are A 2 2 = —A* x 
and A 2 i = A* l2 . These relations imply that a column of Q can be reconstructed using 
just two conjugate gradient inversions of A[&). 

We first examined the timeslice propagator Q{t) = 0; x, t) and empirically 

found the following features: 

• For t ^ 0, Re(N u (t)) w Re(Nu(L t — t)), ie. the antifermion propagator is 
related to that of the fermion by time reversal. 

• lm(N n ) w Im(A^ii) w (iv r 12 ) ~ (iV 12 ) ~ 0: the vanishing of the off-diagonal 
components of N is consistent with isopin SU(2)y symmetry. 
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• Im(A 12 (t)) ~ Im(A 12 (t)) « — Im(A 12 (Li — £)), ie. in the anomalous sector 
fermion and antifermion have equivalent behaviour under time-reversal. 

• Ke(Ai2) ~ (^4ii) ~ (Ali) ~ 0; i e - isopin symmetry in the anomalous sector 
demands the diagonal components vanish. 

From here on for convenience we will denote ReiVn by iV and Iim4 12 by A. In 
Fig. |T^| we plot In |iV(£)| vs. t for two values of \x in the low density chirally broken 
phase. At \x = the propagator is symmetric under time reversal, and the fermion 
mass X ~ 0.730(2), consistent with the breaking of chiral symmetry. At fi — 0.2 the 
time-reversal symmetry is broken; one state propagates forwards with mass 0.952(4), 
the other backwards with mass 0.530(2), corresponding approximately to masses E± 
/i. Now, since tib = for // < /i c , the ground state is unchanged and the physical 
interpretation of this result is simply that the chemical potential shifts the energies 
required to excite fermions and anti-fermions in opposite directions, the anti-fermion 
travelling in the +t direction and the particle —t. 

For n > ji c the situation is completely altered. Recall that in the presence of 
a Fermi surface excitations have a characteristic momentum scale kp. Therefore it 
is necessary to introduce momentum dependence into the Gor'kov propagator via 
Q(k,t) = J2xG{0,0; x,t)e~ zk ^ [jTE| . We choose k oriented along a lattice axis, and 
the set x to include only sites an even number of lattice spacings from the origin in 
each direction, so that the physically accessible momenta are given by k = ^ with 
n = 0, 1, . . . , L s /4. Fig. [13] shows both normal and anomalous propagators in the 
high density phase for k = j. Note that now N(t) ~ for t even, and A(t) ~ for t 
odd. This is a manifestation of the restored SU(2) £ <g> SU(2) B symmetry ( |2.5|) , which 
would be broken by any N e£y00 or A oe eo ^ 0. We have found that the propagators for 
various j and k are well fitted on every second timeslice by the following forms, with 
fit parameters A, B, C and E: 

N(k, t) = Ae~ Et + Be- E{Lt - l \ (5.2) 
A(k,t) = C{e~ Et -e- E{Lt - l) ). (5.3) 

The resulting E(k,j) is the quasiparticle dispersion relation. We have studied this 
function in detail on a 32 3 lattice, which has 9 distinct values of k, performing fits 
to (|5.2j ) over the range t G [5,27]. The fits are of excellent quality, with x 2 /dof 
rarely exceeding 2.0. We are naturally interested in the U(1)b symmetric limit; 



Fig. [14] shows that E(k,j) can be smoothly extrapolated to j = 0, and we quote the 
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Figure 12: Normal propagators |iV(t)| for /i = 0.0 and \i = 0.2 on a 16 2 x 24 lattice 
with j = 0.1. 
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Figure 13: Normal and anomalous propagators with k = ? at = 0.8 on a 16 2 x 24 
lattice with j = 0.1. 
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Figure 14: E(k,j) vs. j extracted from fits to ( |5.2|) on a 32 3 lattice with \x = 0.8. 
Open symbols denote excitations on the hole branch, and closed symbols the particle 
branch. 



results of linearly extrapolating the data with j e [0.025, 0.1]. Fits to the anomalous 
propagator A(k,t) yield quantitatively very similar results for E(k,j). 

The caption of Fig. [14] assigns different k ranges to a "hole branch" (E decreases 
with k) and a "particle branch" (E increases with k). It is straightforward to verify 
this interpretation by considering the free Euclidean propagator Sp(k) = (ift +/i7o + 



m) 1 . For fixed spatial momentum, with jj, < E{k) = yk 2 + m 2 , 

( JZLfl + y-\ e -{E+n)t + > n- 
SAk,t) = {^ + J r >l_ (B _„ M (5.4) 

where the complex "4- velocity" = (E, ±ik)/m. The propagator has both forwards 
and backwards-decaying signals, each associated with a different projection operator; 
in the limit k — > these become ~(1 ± 70) and thus project onto anti-fermion and 
fermion states respectively. The fermion, being lighter, dominates the signal for 
fit ^> 1 yielding a predominantly backward propagation. For fi > E(k), however, 
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there is only a forwards-moving signal, once again dominated by the fermion: 

S F (k, t) = ^Q(t) [(1 + f -) e -^* - (1 - f+) e -^-m . (5.5) 



For free fermions the transition between (|5.4j ) and (|5.5|) takes place at a sharply 



defined Fermi energy Ep = ykp + m 2 = \i. Excitations with k > kp, the Fermi 
momentum, add particles to levels above the Fermi surface; those with k < kp vacate 
holes in the Fermi Sea. The energy cost is smallest when k ~ kp. 

In Fig. [15] we plot the amplitudes A, B and C from the fits to G5.2|J).3|) and 



confirm that for small k, N(k, t) is dominated by a forwards-moving signal, but 
there is a rather sharp crossover to backwards propagation at k/ir ~ 0.3. This 
transition becomes sharper as j — > 0; however, we plot data with j ' ^ to show 
that the amplitude C only differs significantly from zero for momentum states in 
the neighbourhood of the Fermi surface. Were a BCS gap to form, we would expect 
lim^o C(j) indicating particle-hole mixing; our data, however, do not give strong 
support for this. 



Fig. shows the dispersion relation E(k) extrapolated to j = for /i = 0.8, 
together with points derived from the free Gor'kov propagator (ie. generated with 
g 2 = 0) using the identical procedure. We have plotted energies from the hole branch 
as negative in order to generate a smooth curve. There is no sign of any discontinuity 
characteristic of a BCS gap A / 0. In order to interpret the detailed form of the 
curve it is necessary to take account of discretisation effects; for free massless fermions 
the expected dispersion relation, shown as a solid curve in Fig. [16| is E{k) = —fi + 
sinh _1 (sin k). We have found that a reasonable fit to our data for ji 6 [0.8, 0.9], shown 



as a dashed line in Fig. [16], is given by 

E(k) = -E + D smb.' 1 (sink). (5.6) 
Eqn. ( |5.6| ) predicts a sharply defined effective Fermi momentum given by 

K F = sinh _1 (sin k F ) = E /D. (5.7) 

In addition it is possible to define a quasiparticle group velocity (3 = (9sinh(£' + 
E )/dsink, whose value 

coshffp 

Pf = D — — — 5.8 
cosh K p 

at the Fermi surface is the Fermi velocity, which helps to characterise the Fermi 
liquid. For free massless fermions Kp = /i and Pf = 1 for all fi. Our fitted values of 
Kp and Pf are given in Table [2|. 
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Figure 15: The amplitudes A, B and C extracted from fits to ( |5.2| , [5.3|) on a 32 3 
lattice with /i = 0.8. The data for A and B were taken with j = 0.025. 
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Figure 16: Dispersion relation E(k) at fi = 0.8 on a 32 3 lattice for both interacting 
and free fermions. 
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Table 2: Quasiparticle parameters resulting from fits of (|5.6|) to data from a 32 3 
lattice. The quoted errors are purely statistical. 





K F 


Pf 


K F /ii ,p F 


0.80 


0.720(3) 


0.670(3) 


1.34(1) 


0.82 


0.738(3) 


0.671(3) 


1.34(1) 


0.84 


0.773(3) 


0.686(3) 


1.34(1) 


0.86 


0.791(5) 


0.673(4) 


1.37(1) 


0.88 


0.811(5) 


0.628(4) 


1.47(1) 


0.90 


0.836(4) 


0.704(4) 


1.32(1) 



Although the errors quoted in Table |2] are almost certainly underestimated, some 
systematic features are apparent. The observed values of Kp increase smoothly with 
H, and are ~ 90% of their free-field values. The Fermi velocity /3p, by contrast, is more 
or less independent of fi, and only ps 70% of the free-field value; ie. the quasiparticles 
travel at less than the speed of light. In non-relativistic Fermi liquid theory |^T], P2"| , 
the ratio Kp/f3p defines a quantity called the effective mass M*, which need not 
coincide with the mass of the fundamental atomic species M; eg. for the archetypal 
Fermi liquid 3 He in the sub-kelvin (but non-superfluid) regime, M* ~ 3M [p2][ . In 
Landau's theory the ratio for a two-dimensional fluid is given in terms of the dipole 
component of the interaction between quasiparticles, ie. 



— = l + Np cos 0, (5.9) 

M Jo Ztt 



where Np = qVv-^S^ is the number of quasiparticle states on the Fermi surface per 
unit energy interval (V-j is the volume of 2c? space and g counts independent spin 
and isospin components) and /($) the spin-singlet interaction energy between quasi- 
particles at the Fermi surface with momenta separated by angle In a relativistic 
generalisation the left hand side of (|5.9|) is replaced by Kp/fi/3p p8f; we infer from 



the data of Table |2] that /($) costf > when averaged over the Fermi circle, and the 
interaction hence repulsive between quasiparticles with parallel momenta and/or at- 
tractive if the momenta are anti-parallel (the simpler conclusion that the interaction 
is always attractive was wrongly drawn in |]I6| ). This should be contrasted with the 
interaction between the fundamental fermions of the NJL model due to a exchange, 
which is attractive and independent of direction^. A similar effective reversal of sign 



arises in the Hartree-Fock treatment of free electron states in a metal [131, and is 



1 Single 7r exchange is attractive in isosinglet but repulsive in isotriplet channels; the net binding 
effect in matter made from equal numbers of u and d quarks vanishes. 



32 



characteristic of a quantum-mechanical exchange effect. 

To summarise, we have examined the quasiparticle spectrum and estimated both 
Fermi momentum Kp and velocity (5p with due allowance made for discretisation 
artifacts. The results are consistent with a relativistic generalisation of a Landau 
Fermi liquid, and qualitatively similar to the normal phase of liquid 3 He. There is 
no evidence for a BCS gap, and in the j — > limit the anomalous components of 
the propagator signalling particle-hole mixing probably vanish. We note, however, 
that superfluid behaviour is not precluded by the absence of a gap [^2|; the unlimited 
growth of quasiparticle excitations that couple normal and superfluid components 
and hence destroy superfluidity in a gapless Bose liquid are here prevented by the 
Pauli Exclusion Principle. Long range phase coherence is a sufficient condition for 
superfluidity. 

6 Summary and Outlook 

Let us briefly review the main achievements of the paper. Firstly, we have developed 
the necessary formalism to identify diquark condensation in numerical lattice studies 
of field theories on finite systems at non-zero chemical potential, the crucial ingredient 
being the introduction of a diquark source term. Secondly, to our initial surprise, 
we have found no evidence for a condensate (qq) ^ in studies of the 2+ld NJL 
model in its high density phase \x > \x c . Rather, the results from two independent 
analyses of diquark observables are consistent with a critical behaviour (qq(j)) oc j a 
throughout the dense phase. Whilst there is some residual uncertainty about the 
source of finite volume effects for j < 0.05, we suspect it would require computer 
resources considerably greater than those we have used to modify this conclusion. 
Critical behaviour in two dimensions implies long range coherence in the phase of the 
condensate wavefunction, which is a sufficient condition for superfluidity. Whilst there 
is qualitative similarity with the well-known example of the low temperature phase of 
the 2d XY model, the measured value of the critical exponent 5 is distinct, suggesting 
that the 2+ld dense NJL model belongs to a new universality class. Thirdly, we have 
performed the first systematic spectroscopic study in the spin-^ sector for \i ^ 0, and 
mapped out the quasiparticle dispersion relation. The success of the simple pole fits 
( |5.2|, |5.3|) confirms the long-lived nature of the quasiparticles. There is no evidence for 
either particle- hole mixing or A ^ in the j — > limit. Instead, the system resembles 
a normal Fermi liquid with a well-defined Fermi surface; the Fermi velocity f3p is of 
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the same order as but significantly less than the free-field value 1. Our findings can 
be summarised by the statement that the high density phase of the 2+ld NJL model 
is a relativistic gapless thin film BCS superfluid. 

In a sense this and related papers with /i ^ represent the primitive beginnings 
of the study of condensed matter physics on the lattice. Let us sketch a few possible 
future directions. Firstly it would be interesting to estimate the supercurrent J s cor- 
responding to a quantised flow pattern around a finite system as in ([4.5|), which could 
be set up using a spatially varying j(x). As well as providing a direct demonstration 
of superfluidity, this would also enable the extraction of the phenomenologically im- 
portant parameter K s . Secondly, it is possible to study quasiparticles and other Fermi 
surface-related phenomena. For instance, a sharp Fermi surface leads to oscillations 
of spatial frequency 2k p in the screened potential between static charges, known as 



Friedel oscillations |13|, |39|. Friedel oscillations can be observed in the wavefunctions 
of qq and qq states in 2+ld four-fermi models with /x > \i c fllOfl . Another possibility is 
the observation of light qq mesons in the spin-1 channel, corresponding to low-energy 



excitations of the Fermi surface related to the phenomenon of zero sound |Ti], |22|, f4T |. 
Finally, it is of prime importance to extend our calculations to the physically relevant 
case of 3+ Id. In this case the NJL model is no longer a fundamental field theory, but 
instead can be thought of as an effective description of strong interaction physics with 
many possible phenomenological applications including thermodynamics |42| . If our 
arguments are correct, then the "conventional" signals for a BCS mechanism, namely 
(qq) 7^ and A^O should be readily observed using the methods we have developed. 
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